Análisis geospacial con Python

Introducción

Paquetes de análisis espacial:

• Geopandas: paquete para datos vectoriales (equivalente a sf)

• Xarray: paquete para datos multidimensionales (equivalente a stars?)

• Rasterio: interfaz a GDAL para datos raster (equivalente a terra, pero de bajo nivel)

• Rioxarray: puente entre Xarray y Rasterio.

Paquetes de visualización interactiva:

• Leafmap: interfaz sencilla y unificada para mapas web de múltiples backends (folium, ipyleaflet, maplibre).

• MapLibre: mapas 2D y 3D avanzados.

• HyperCoast: visualización 3D de datos hiperespectrales con pocas líneas de código.

Paquetes de análisis especializado:

• WhiteboxTools: herramientas de geoprocesamiento de análisis hidrológico, análisis del terreno, y análisis de datos LiDAR.

• Geemap: puente a Google Earth Engine.

• DuckDB: análisis de datos espaciales en base de datos de alto rendimiento.

• Apache Sedona: computación distribuida para procesado de datos geoespaciales.

Desarrollo de aplicaciones:

• Voila y Solara: creación de aplicaciones con Python y Jupyter.

• Shiny: creación de aplicaciones y dashboards.

Vamos a comprobar que los paquetes se cargan correctamente y a crear un primer mapa:

## cargar paquetes
import geopandas as gpd
import rasterio
import xarray as xr
import rioxarray
import leafmap
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

## crear mapa
map = leafmap.Map(center = [40, -100], zoom = 5, height = '500px')

## añadir mapas base
map.add_basemap('OpenTopoMap')
map.add_basemap('USGS.Imagery')

## mostrar mapa
map

Análisis vectorial con GeoPandas

Vamos aprender:

• GeoDataFrame y GeoSeries

• Leer/Exportar datos vectoriales

• Operaciones comunes

• Visualizar datos

• Trabajar con CRS

Geopandas

Una diferencia con sf, es que un GeoDataFrame puede tener varias geometrías, pero solamente una está activa en cada momento. La geometría se accede a través del atributo .geometry. Vamos a crear un GeoDataFrame:

## datos
data = {
    'City': ['Tokyo', 'New York', 'London', 'Paris'],
    "Latitude": [35.6895, 40.7128, 51.5074, 48.8566],
    "Longitude": [139.6917, -74.0060, -0.1278, 2.3522]
}

## convertir a DataFrame
df = pd.DataFrame(data)

## convertir a GeoDataFrame
gdf = gpd.GeoDataFrame(df, geometry = gpd.points_from_xy(df.Longitude, df.Latitude))

## imprimir
print(gdf)

       City  Latitude  Longitude                  geometry
0     Tokyo   35.6895   139.6917  POINT (139.6917 35.6895)
1  New York   40.7128   -74.0060   POINT (-74.006 40.7128)
2    London   51.5074    -0.1278   POINT (-0.1278 51.5074)
3     Paris   48.8566     2.3522    POINT (2.3522 48.8566)

Leer datos:

## url a geojson
url = "https://github.com/opengeos/datasets/releases/download/vector/nybb.geojson"

## leer
data_geojson = gpd.read_file(url)

Exportar datos:

gdf.to_file('00_data\\test-data.geojson')
gdf.to_file('00_data\\test-data.gpkg', layer = "layer1")

gpd.list_layers('00_data\\test-data.gpkg')

	name	geometry_type
0	layer1	Point

Sistemas de Referencia de Coordenadas

Comprobar CRS:

print(data_geojson.crs)

EPSG:2263

Reproyectar o transformar CRS:

gdf_4326 = data_geojson.to_crs('EPSG:4326')

Medidas

Para calular medidas como áreas y distancias, necesitamos un CRS que esté en unidades como metros (p. ej. EPSG:3857).

## reproyectar
gdf_3857 = gdf_4326.to_crs('EPSG:3857')

## crear índice
gdf_3857.set_index('BoroName', inplace = True)

Calcular área:

## área
gdf_3857["area"] = gdf_3857.area

## perímetro en km
gdf_3857["perim"] = gdf_3857.length / 1000

Extraer características geométricas

## centroide
gdf_3857['centroid'] = gdf_3857.centroid

## bordes
gdf_3857['boundary'] = gdf_3857.boundary

Distancias: vamos a ver la distancia desde el centroide de Manhattan al resto de centroides

## centroide de referencia
manhattan_centroid = gdf_3857.loc['Manhattan', 'centroid']

## calcular distancia
gdf_3857['dist_to_manhattan'] = gdf_3857['centroid'].distance(manhattan_centroid)

Crear buffers:

## crear buffer (metros)
gdf_3857['buffer'] = gdf_3857.buffer(3000)

## visualizar
fig, ax = plt.subplots(figsize=(10, 6))

gdf_3857["buffer"].plot(
    ax        = ax,
    alpha     = 0.3,
    color     = "orange",
    edgecolor = "red",
    linewidth = 1,
    label     = "3km Buffer Zone",
)

gdf_3857.plot(
    ax        = ax,
    color     = 'lightblue',
    linewidth = 1,
    edgecolor = 'navy',
    label     = 'Original Boundaries'
)

plt.title('Distritos y Buffer')
plt.legend()
plt.axis('off')

(np.float64(-8272486.092246463),
 np.float64(-8197855.319316324),
 np.float64(4931921.070848358),
 np.float64(5006269.877978747))

Crear convex hulls:

gdf_3857['chull'] = gdf_3857.convex_hull

Visualización

Vamos a especificar parámetros para mejorar las visualizaciones:

plt.rcParams['figure.dpi'] = 150

Mapa temático:

## crear lienzo
fig, ax = plt.subplots(figsize = (10, 6))

## añadir mapa
gdf_3857.plot(
    column    = 'area',
    ax        = ax,
    legend    = True,
    cmap      = 'YlOrRd',
    edgecolor = 'black',
    linewidth = .5
)

## añadir elementos al mapa
plt.title('Distritos de NYC por área', fontsize = 16, fontweight = 'bold')
plt.axis('off')
plt.tight_layout()
plt.show()

Mapa interactivo con Folium:

gdf_3857.explore(
    column  = 'area',
    cmap    = 'YlOrRd',
    tooltip = ['area', 'dist_to_manhattan'],
    popup   = True,
    legend  = True
)

Make this Notebook Trusted to load map: File -> Trust Notebook

Relaciones espaciales

Las relaciones espaciales en Python se estrucutran: objeto1.relación(objeto2). Los métodos intersects, touches, crosses, contains, contains_properly… devuelven un Pandas Series de True/False.

Intersección con Manhattan:

## extraer Manhattan
manhattan_geom = gdf_3857.loc['Manhattan', 'geometry']

## intersección
intersects_manhattan = gdf_3857.intersects(manhattan_geom)

## filtrar los que intersecan
gdf_3857.loc[intersects_manhattan, :]

	BoroCode	Shape_Leng	Shape_Area	geometry	area	perim	centroid	boundary	dist_to_manhattan	buffer	chull
BoroName
Queens	4	896344.047763	3.045213e+09	MULTIPOLYGON (((-8219461.955 4952778.645, -821...	4.928308e+08	360.377573	POINT (-8217436.761 4969318.623)	MULTILINESTRING ((-8219461.955 4952778.645, -8...	19456.427458	POLYGON ((-8229469.036 4965425.107, -8229556.9...	POLYGON ((-8231045.21 4944993.298, -8233479.86...
Brooklyn	3	741080.523166	1.937479e+09	MULTIPOLYGON (((-8222843.701 4950893.717, -822...	3.129684e+08	297.784594	POINT (-8231817.476 4960085.209)	MULTILINESTRING ((-8222843.701 4950893.717, -8...	19586.763567	POLYGON ((-8245284.602 4956687.136, -8245296.1...	POLYGON ((-8235790.052 4949053.259, -8237874.4...
Manhattan	1	359299.096471	6.364715e+08	MULTIPOLYGON (((-8238858.852 4965914.967, -823...	1.032201e+08	144.706085	POINT (-8233984.776 4979551.696)	MULTILINESTRING ((-8238858.852 4965914.967, -8...	0.000000	POLYGON ((-8237068.209 4963506.169, -8237106.8...	POLYGON ((-8240208.861 4965683.834, -8242627.6...
Bronx	2	464392.991824	1.186925e+09	MULTIPOLYGON (((-8226155.114 4982269.852, -822...	1.929250e+08	187.221157	POINT (-8222783.625 4990631.161)	MULTILINESTRING ((-8226155.114 4982269.852, -8...	15755.010025	POLYGON ((-8233207.11 4988220.645, -8233189.61...	POLYGON ((-8224095.48 4980733.074, -8225054.51...

Tocan a Manhattan:

gdf_3857.loc[gdf_3857.touches(manhattan_geom), :]

	BoroCode	Shape_Leng	Shape_Area	geometry	area	perim	centroid	boundary	dist_to_manhattan	buffer	chull
BoroName
Queens	4	896344.047763	3.045213e+09	MULTIPOLYGON (((-8219461.955 4952778.645, -821...	4.928308e+08	360.377573	POINT (-8217436.761 4969318.623)	MULTILINESTRING ((-8219461.955 4952778.645, -8...	19456.427458	POLYGON ((-8229469.036 4965425.107, -8229556.9...	POLYGON ((-8231045.21 4944993.298, -8233479.86...
Brooklyn	3	741080.523166	1.937479e+09	MULTIPOLYGON (((-8222843.701 4950893.717, -822...	3.129684e+08	297.784594	POINT (-8231817.476 4960085.209)	MULTILINESTRING ((-8222843.701 4950893.717, -8...	19586.763567	POLYGON ((-8245284.602 4956687.136, -8245296.1...	POLYGON ((-8235790.052 4949053.259, -8237874.4...
Bronx	2	464392.991824	1.186925e+09	MULTIPOLYGON (((-8226155.114 4982269.852, -822...	1.929250e+08	187.221157	POINT (-8222783.625 4990631.161)	MULTILINESTRING ((-8226155.114 4982269.852, -8...	15755.010025	POLYGON ((-8233207.11 4988220.645, -8233189.61...	POLYGON ((-8224095.48 4980733.074, -8225054.51...

Comprobar que los centroides están dentro de su polígono:

gdf_3857['centroid'].within(gdf_3857.geometry)

BoroName
Staten Island    True
Queens           True
Brooklyn         True
Manhattan        True
Bronx            True
dtype: bool

Relaciones binarias (dos geometrías)

gdf_3857['intersection'] = gdf_3857.intersection(manhattan_geom)
gdf_3857['difference'] = gdf_3857.difference(manhattan_geom)

gdf_3857['intersection'].explore()

Make this Notebook Trusted to load map: File -> Trust Notebook

Verificaciones

gdf.is_empty
gdf.is_valid
gdf.is_simple

0    True
1    True
2    True
3    True
dtype: bool